Package deSolve: Solving Initial Value Differential Equations in R

Demonstration of the package with plotting in ggplot2 and plotly

Introduction

This notebook is a demonstration of the R package deSolve. It is a direct implementation of the code in this vignette provided by the package authors, (Soetaert, Petzoldt and Setzer). However the plots are implemented using ggplot2, with plotly used for some 3-d examples. The plotting functions included in the package use base R. The purpose of this notebook is to demonstrate the use of modern plotting libraries with the output of the deSolve package.

Setup

library(deSolve)
library(tidyverse)
theme_set(theme_bw())
library(patchwork)
library(plotly)
library(metR)

1 A simple ODE: chaos in the atmosphere

We use the Lorenz model to demonstrate how to implement and solve differential equations in R. The Lorenz model describes the dynamics of three state variables, \(X\), \(Y\) and \(Z\). The model equations are: \[ \frac{dX}{dt} = a \cdot X + Y \cdot Z \] \[ \frac{dY}{dt} = b \cdot (Y - Z) \] \[ \frac{dZ}{dt} = -X \cdot Y + c \cdot Y - Z \] with the initial conditions: \[ X(0) = Y(0) = Z(0) = 1 \] where \(a\), \(b\) and \(c\) are three parameters with values of -8/3, -10 and 28 respectively.

Implementation of an IVP ODE in R can be separated in two parts: the model specification and the model application. Model specification consists of:

• Defining model parameters and their values,
• Defining model state variables and their initial conditions,
• Implementing the model equations that calculate the rate of change (e.g. \(dX/dt\)) of the state variables.

The model application consists of:

• Specification of the time at which model output is wanted,
• Integration of the model equations (uses R functions from deSolve),
• Plotting of model results. Below, we discuss the R code for the Lorenz model.

1.1 Model specification

Model parameters

parameters <- c(a = -8/3, b = -10, c = 28)

State variables

state <- c(X = 1, Y = 1, Z = 1)

Model equations

Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    # rate of change
    dX <- a*X + Y*Z
    dY <- b * (Y-Z)
    dZ <- -X*Y + c*Y - Z
  
    # return the rate of change
    list(c(dX, dY, dZ))
  })
}

1.2 Model application

Time specification

times <- seq(0, 100, by = 0.01)

Model integration

out <- ode(y = state,
           times = times,
           func = Lorenz,
           parms = parameters
       )
head(out)

     time         X        Y        Z
[1,] 0.00 1.0000000 1.000000 1.000000
[2,] 0.01 0.9848912 1.012567 1.259918
[3,] 0.02 0.9731148 1.048823 1.523999
[4,] 0.03 0.9651593 1.107207 1.798314
[5,] 0.04 0.9617377 1.186866 2.088545
[6,] 0.05 0.9638068 1.287555 2.400161

Plotting results

Plot the results using ggplot2:

# Convert data to long format:
out_long <- data.frame(out) %>% 
  pivot_longer(cols = X:Z, names_to = "var", values_to = "value")

# Plot X, Y and Z versus time:
P1 <- out_long %>% 
  split(.$var) %>% 
  map(function(df) {
    ggplot(df) +
      geom_line(aes(time, value), color = "steelblue") +
      scale_x_continuous(breaks = seq(0, 100, 20)) +
      labs(x = "time", y = df$var[1],
           title = paste0("<span style=\"font-family:'Lucida Console', monospace\">", df$var[1], "</span> versus <span style=\"font-family:'Lucida Console', monospace\">time</span>")) + # Ref: https://wilkelab.org/ggtext/
      theme(panel.grid = element_blank(),
            plot.title = ggtext::element_markdown())
  })

# Plot Z versus X:
P1[[4]] <- ggplot(data.frame(out)) +
  geom_point(aes(X, Z), shape = '.', color = "steelblue") +
  ggtitle("Z versus X") +
  theme(panel.grid = element_blank())

# Output plots as 2x2 grid:
(P1[[1]] + P1[[2]]) / (P1[[3]] + P1[[4]]) +
  plot_annotation(title = "Figure 1: Solution of the ordinary differential equation")

2 Solvers for initial value problems of ordinary differential equations

Demonstrate the run times of a few different solvers from the deSolve package:

# List of integration routines from deSolve to include:
routines <-
  list(
    rk4 = rk4,
    lsode = lsode,
    lsoda = lsoda,
    lsodes = lsodes,
    daspk = daspk,
    vode = vode
)

# Run each routine and display the run time as a data frame:
args <- list(y = state, times = times, func = Lorenz, parms = parameters)
times_list <- map(routines, function(f) system.time(exec(f, !!!args))) # https://stackoverflow.com/a/61132298
map(times_list, as.vector) %>% 
  dplyr::bind_rows() %>% 
  head(3) %>% 
  cbind(var = c("user", "system", "elapsed")) %>% 
  pivot_longer(cols = !contains("var"), names_to = "routine") %>% 
  pivot_wider(id_cols = routine, names_from = var)

2.1. Runge-Kutta methods and Euler

Available methods:

rkMethod()

 [1] "euler"   "rk2"     "rk4"     "rk23"    "rk23bs"  "rk34f"   "rk45f"   "rk45ck"  "rk45e"   "rk45dp6" "rk45dp7" "rk78dp" 
[13] "rk78f"   "irk3r"   "irk5r"   "irk4hh"  "irk6kb"  "irk4l"   "irk6l"   "ode23"   "ode45"  

Define and use a new Runge-Kutta method, rKnew:

func <- function(t, x, parms) {
  with(as.list(c(parms, x)),{
    dP <- a * P     - b * C * P
    dC <- b * P * C - c * C
    res <- c(dP, dC)
    list(res)
  })
}

rKnew <- rkMethod(ID = "midpoint",
  varstep = FALSE,
  A = c(0, 1/2),
  b1 = c(0, 1),
  c = c(0, 1/2),
  stage = 2,
  Qerr = 1
)

out <- ode(y = c(P = 2, C = 1),
           times = 0:100, func,
           parms = c(a = 0.1, b = 0.1, c = 0.1),
           method = rKnew
       )
head(out)

     time        P        C
[1,]    0 2.000000 1.000000
[2,]    1 1.990000 1.105000
[3,]    2 1.958387 1.218598
[4,]    3 1.904734 1.338250
[5,]    4 1.830060 1.460298
[6,]    5 1.736925 1.580136

2.2. Model diagnostics and summaries

The function diagnostics prints several diagnostics of the simulation to the screen:

out1 <- rk4(state, times, Lorenz, parameters)
out2 <- lsode(state, times, Lorenz, parameters)

diagnostics(out1)

--------------------
rk return code
--------------------

  return code (idid) =  0 
  Integration was successful.

--------------------
INTEGER values
--------------------

  1 The return code : 0 
  2 The number of steps taken for the problem so far: 10000 
  3 The number of function evaluations for the problem so far: 40001 
 18 The order (or maximum order) of the method: 4 
 

diagnostics(out2)

--------------------
lsode return code
--------------------

  return code (idid) =  2 
  Integration was successful.

--------------------
INTEGER values
--------------------

  1 The return code : 2 
  2 The number of steps taken for the problem so far: 12778 
  3 The number of function evaluations for the problem so far: 16633 
  5 The method order last used (successfully): 5 
  6 The order of the method to be attempted on the next step: 5 
  7 If return flag =-4,-5: the largest component in error vector 0 
  8 The length of the real work array actually required: 58 
  9 The length of the integer work array actually required: 23 
 14 The number of Jacobian evaluations and LU decompositions so far: 721 
 
--------------------
RSTATE values
--------------------

  1 The step size in t last used (successfully): 0.01 
  2 The step size to be attempted on the next step: 0.01 
  3 The current value of the independent variable which the solver has reached: 100.0072 
  4 Tolerance scale factor > 1.0 computed when requesting too much accuracy: 0 
 

There is also a summary method for deSolve objects:

summary(out1)

3 Partial differential equations

Consider the Aphid model described in Soetaert and Herman (2009)¹. It is a model where aphids (a pest insect) slowly diffuse and grow on a row of plants. The model equations are: \[ \frac{\partial N}{\partial t}=-\frac{\partial Flux}{\partial x} + g \cdot N \] and where the diffusive flux is given by: \[ Flux = -D \frac{\partial N}{\partial x} \] with boundary conditions \[ N_{x=0} = N_{x=60} = 0 \] and initial condition \[ N_x = 0 \text{ for } x \ne 30 \] \[ N_x = 1 \text{ for } x = 30 \] In the method of lines approach, the spatial domain is subdivided in a number of boxes and the equation is discretized as: \[ \frac{dN_i}{dt} = - \frac{Flux_{i,i+1} - Flux_{i-1,i}}{\Delta x_i} + g \cdot N_i \] with the flux on the interface equal to: \[ Flux_{i-1,i} = -D_{i-1,i} \cdot \frac{N_i - N_{i-1}}{\Delta x_{i-1,i}} \]

Define the model equations:

Aphid <- function(t, APHIDS, parameters) {
  deltax <- c(0.5, rep(1, numboxes - 1), 0.5)
  Flux <- -D * diff(c(0, APHIDS, 0)) / deltax
  dAPHIDS <- -diff(Flux) / delx + APHIDS * r
  
  # Return:
  list(dAPHIDS)
}

Model parameters:

D <- 0.3    # m2/day   diffusion rate
r <- 0.01   # /day     net growth rate
delx <- 1   # m        thickness of boxes
numboxes <- 60

# distance of boxes on plant, m, 1m intervals
Distance <- seq(from = 0.5, by = delx, length.out = numboxes)

Aphids are initially only present in two central boxes:

# Initial condistions:   # ind/m2
APHIDS <- rep(0, times = numboxes)
APHIDS[30:31] <- 1
state <- c(APHIDS = APHIDS) # initialise state variables

The model is run for 200 days, producing output every day; the time elapsed in seconds to solve this 60 state-variable model is estimated with system.time:

times <- seq(0, 200, by = 1)
print(system.time(
  out <- ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid")
))

   user  system elapsed 
  0.039   0.000   0.040 

The matrix out consist of times (1st column) followed by the densities (next columns):

head(out[ , 1:5])

     time      APHIDS1      APHIDS2      APHIDS3      APHIDS4
[1,]    0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[2,]    1 1.667194e-55 9.555028e-52 2.555091e-48 4.943131e-45
[3,]    2 3.630860e-41 4.865105e-39 5.394287e-37 5.053775e-35
[4,]    3 2.051210e-34 9.207997e-33 3.722714e-31 1.390691e-29
[5,]    4 1.307456e-30 3.718598e-29 9.635350e-28 2.360716e-26
[6,]    5 6.839152e-28 1.465288e-26 2.860056e-25 5.334391e-24

Summary of the entire 1-D variable:

summary(out)

Plot the output using ggplot2

First, a contour plot:

# Create a data frame of APHID numbers and distances for joining to the output:
Distance_APHIDS <- data.frame(
  APHIDS_ = paste0("APHIDS", 1:length(Distance)),
  Distance
)

# Pivot the output to longer format and join with the above data frame to add the Distance values:
out_long <- data.frame(out) %>% 
  pivot_longer(cols = starts_with("APHIDS"), names_to = "APHIDS_", values_to = "density") %>%
  left_join(Distance_APHIDS, by = "APHIDS_")

# Plot:
ggplot(out_long) +
  geom_raster(aes(time, Distance, fill = density),
             interpolate = TRUE) +
  geom_contour(aes(time, Distance, z = density), bins = 12,
               color = "white", size = 0.1, alpha = 0.5) +
  scale_fill_distiller(palette = "Spectral",
                       breaks = seq(0, 1, 0.2)) +
  labs(x = "Time (days)", y = "Distance on plant (m)", fill = "Aphid density",
       title = "Figure 2: Solution of the 1-dimensional aphid model",
       subtitle = "Aphid density on a row of plants") +
  guides(fill = guide_colourbar(barheight = unit(3.3, "in"),
                                ticks.colour = "black",
                                frame.colour = "black")) +
  coord_cartesian(expand = F)

Line plots:

points <-
  data.frame(
    Distance = c(0, 10, 20, 30, 40, 50, 60),
    density = c(0, 0.1, 0.25, 0.5, 0.25, 0.1, 0)
  )

p1 <- out_long %>% 
  filter(time %in% seq(0, 200, by = 10)) %>% 
  ggplot() +
  geom_line(aes(Distance, density, group = time, linetype = factor(time)),
            color = "steelblue", show.legend = FALSE) +
  geom_point(data = points, aes(Distance, density),
             color = "indianred", size = 3) +
  scale_x_continuous(breaks = seq(0, 60, 10)) +
  scale_y_continuous(breaks = seq(0, 1, 0.2)) +
  labs(x = "x", y = "", title = "Density at 21 evenly space times from 0 to 200") +
  theme(panel.grid = element_blank(),
        plot.title = element_text(size = 12))

p2 <- out_long %>% 
  filter(time == 100) %>% 
  ggplot() +
  geom_line(aes(Distance, density, group = time, linetype = factor(time)),
            color = "steelblue", show.legend = FALSE) +
  geom_point(data = points, aes(Distance, density),
             color = "indianred", size = 3) +
  scale_x_continuous(breaks = seq(0, 60, 10)) +
  scale_y_continuous(breaks = seq(0, 0.5, 0.1)) +
  labs(x = "x", y = "", title = paste0("<span>Density at </span><span style=\"font-family:'Lucida Console', monospace\">time = 100</span>")) +
  theme(panel.grid = element_blank(),
        plot.title = ggtext::element_markdown(size = 12))

p1 + p2 + plot_annotation(title = "Figure 3: Solution of the Aphid model")

4 Differential algebraic equations

4.1. DAEs of index maximal 1

Consider the following simple DAE: \[ \frac{dy_1}{dt} = -y_1 + y_2 \] \[ y_1 \cdot y_2 = t \] where the first equation is a differential, the second an algebraic equation. To solve it, it is first rewritten as residual functions: \[ 0 = \frac{dy_1}{dt} + y_1 - y_2 \] \[ 0 = y_1 \cdot y_2 - t \]

In R we write:

daefun <- function(t, y, dy, parameters) {
  res1 <- dy[1] + y[1] - y[2]
  res2 <- y[2] * y[1] - t
  
  # Return:
  list(c(res1, res2))
}

yini <- c(1, 0)
dyini <- c(1, 0)
times <- seq(0, 10, 0.1)

system.time(out <- daspk(y = yini, dy = dyini, times = times, res = daefun, parms = 0))

   user  system elapsed 
  0.012   0.000   0.012 

Plot using ggplot2:

pendulum <- function (t, Y, parms) {
  with (as.list(Y),
  list(c(u,
         v,
         -lam * x,
         -lam * y - 9.8,
         x^2 + y^2 -1
  ))
 )
}

4.2. DAEs of index up to three

Consider the well-known pendulum equation:

\[ x' = u \] \[ y' = v \] \[ u' = - \lambda x \] \[ v' = -\lambda y - 9.8 \] \[ 0 = x^2 + y^2 - 1 \]

where the dependent variables are \(x\), \(y\), \(u\), \(v\) and \(\lambda\).

Implemented in R to be used with function radau this becomes:

pendulum <- function (t, Y, parms) {
  with (as.list(Y),
  list(c(u,
         v,
         -lam * x,
         -lam * y - 9.8,
         x^2 + y^2 -1
  ))
 )
}

A consistent set of initial conditions are:

yini <- c(x = 1, y = 0, u = 0, v = 1, lam = 1)

Mass matrix:

M <- diag(nrow = 5)
M[5, 5] <- 0
M

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    0

Function radau requires that the index of each equation is specified; there are 2 equations of index 1, two of index 2, one of index 3:

index <- c(2, 2, 1)
times <- seq(from = 0, to = 10, by = 0.01)
out <- radau (y = yini, func = pendulum, parms = NULL, times = times, mass = M, nind = index)

Plot using ggplot2:

out_long <- data.frame(out) %>% 
  pivot_longer(cols = x:lam, names_to = "var")

# Plot X, Y and Z versus time:
P5 <- out_long %>% 
  split(.$var) %>% 
  map(function(df) {
    ggplot(df) +
      geom_line(aes(time, value), color = "steelblue") +
      scale_x_continuous(breaks = seq(0, 100, 20)) +
      labs(x = "time", y = df$var[1],
           title = paste0("<span style=\"font-family:'Lucida Console', monospace\">", df$var[1], "</span> versus <span style=\"font-family:'Lucida Console', monospace\">time</span>")) + # Ref: https://wilkelab.org/ggtext/
      theme(panel.grid = element_blank(),
            plot.title = ggtext::element_markdown(size = 12))
  })

# Plot y versus x:
P5[[6]] <- ggplot(data.frame(out)) +
  geom_line(aes(x, y), color = "steelblue") +
  ggtitle("y versus x") +
  theme(panel.grid = element_blank(),
        plot.title = element_text(size = 12))

# Output plots as 2x3 grid:
(P5[[4]] + P5[[5]] + P5[[2]]) / (P5[[3]] + P5[[1]] + P5[[6]]) +
  plot_annotation(title = "Figure 5: Solution of the pendulum problem, an index 3 differential algebraic equation")

7 Events and roots

7.1 Event specified in a data.frame

In this example, two state variables with constant decay are modeled:

eventmod <- function(t, var, parms) {
  list(dvar = -0.1*var)
}
yini <- c(v1 = 1, v2 = 2)
times <- seq(0, 10, by = 0.1)

At time 1 and 9 a value is added to variable v1, at time 1 state variable v2 is multiplied with 2, while at time 5 the value of v2 is replaced with 3. These events are specified in a data.frame called eventdat:

ballode<- function(t, y, parms) {
  dy1 <- y[2]
  dy2 <- -9.8
  
  # Return:
  list(c(dy1, dy2))
}

The model is solved with ode:

out <- ode(func = eventmod, y = yini, times = times, parms = NULL, events = list(data = eventdat))

out_long <- data.frame(out) %>% 
  pivot_longer(cols = v1:v2, names_to = "var")

P6 <- out_long %>% 
  split(.$var) %>% 
  map(function(df) {
    ggplot(df) +
      geom_line(aes(time, value), color = "steelblue") +
      scale_x_continuous(breaks = seq(0, 10, 2)) +
      labs(x = "time", y = df$var[1],
           title = paste0("<span style=\"font-family:'Lucida Console', monospace\">", df$var[1], "</span> versus <span style=\"font-family:'Lucida Console', monospace\">time</span>")) + # Ref: https://wilkelab.org/ggtext/
      theme(panel.grid = element_blank(),
            plot.title = ggtext::element_markdown(size = 12))
  })

P6[[1]] + P6[[2]]  +
  plot_annotation(title = "Figure 6: A simple model that contains events")

7.2 Event triggered by a root function

This model describes the position (y1) and velocity (y2) of a bouncing ball:

ballode<- function(t, y, parms) {
  dy1 <- y[2]
  dy2 <- -9.8
  
  # Return:
  list(c(dy1, dy2))
}

An event is triggered when the ball hits the ground (height = 0) Then velocity (y2) is reversed and reduced by 10 percent. The root function, y[1] = 0, triggers the event:

root <- function(t, y, parms) y[1]

The event function imposes the bouncing of the ball:

event <- function(t, y, parms) {
  y[1]<- 0
  y[2]<- -0.9 * y[2]
  return(y)
}

After specifying the initial values and times, the model is solved, here using lsode:

yini <- c(height = 0, v = 20)
times <- seq(from = 0, to = 20, by = 0.01)
out <- lsode(times = times, y = yini, func = ballode, parms = NULL,
             events = list(func = event, root = TRUE),
             rootfun = root)

data.frame(out) %>% 
  ggplot() +
  geom_line(aes(time, height), color = "steelblue") +
  ggtitle("Figure 7: A model, with event triggered by a root function") +
  theme(panel.grid = element_blank())

8 Delay differential equations

We implement the lemming model, example 6 from (Shampine and Thompson 2000)². Function lagvalue calculates the value of the state variable at t - 0.74. As long a these lag values are not known, the value 19 is assigned to the state variable. Note that the simulation starts at time = - 0.74:

# the derivative function
derivs <- function(t, y, parms) {
  if (t < 0)
    lag <- 19
  else
    lag <- lagvalue(t - 0.74)
  
  dy <- r * y * (1 - lag/m)
  list(dy, dy = dy)
}

# parameters
r <- 3.5; m <- 19

# initial values and times
yinit <- c(y = 19.001)
times <- seq(-0.74, 40, by = 0.01)

# solve the model
yout <- dede(y = yinit, times = times, func = derivs, parms = NULL, atol = 1e-10)

# plot using ggplot2
P8a <- ggplot(data.frame(yout)) +
  geom_line(aes(time, y), color = "steelblue") +
  scale_y_continuous(breaks = seq(0, 100, 20)) +
  theme(panel.grid = element_blank())

P8b <- ggplot(data.frame(yout)) +
  geom_path(aes(y, dy.y), color = "steelblue") +
  scale_x_continuous(breaks = seq(0, 100, 20)) +
  scale_y_continuous(breaks = seq(-200, 100, 100)) +
  labs(y = "dy") +
  theme(panel.grid = element_blank())

P8a + P8b + plot_annotation(title = "Figure 8: A delay differential equation model (Lemming model)")

9 Discrete time models, difference equations

We give here an example of a discrete time model, represented by a difference equation: the Teasel model as from Soetaert and Herman (2009, p287). The dynamics of this plant is described by 6 stages and the transition from one stage to another is in a transition matrix.

We define the stages and the transition matrix first:

Stages <- c("DS 1yr", "DS 2yr", "R small", "R medium", "R large", "F")
NumStages <- length(Stages)

# Population matrix
A <- matrix(nrow = NumStages, ncol = NumStages, byrow = TRUE,
            data = c(0,     0,    0,     0,     0,     322.38,
                     0.966, 0,    0,     0,     0,     0,
                     0.013, 0.01, 0.125, 0,     0,     3.448,
                     0.007, 0,    0.125, 0.238, 0,     30.170,
                     0.008, 0,    0.038, 0.245, 0.167, 0.862,
                     0,     0,    0,     0.023, 0.75,  0)
           )

The difference function is defined as usual, but does not return the “rate of change” but rather the new relative stage densities are returned. Thus, each time step, the updated values are divided by the summed densities:

Teasel <- function (t, y, p) {
  yNew <- A %*% y
  list (yNew / sum(yNew))
}

The model is solved using method "iteration":

out <- ode(func = Teasel, y = c(1, rep(0, 5) ), times = 0:50, parms = 0, method = "iteration")

Plot using ggplot2:

data.frame(out) %>% 
  pivot_longer(cols = X1:X6, names_to = "var") %>% 
  ggplot() +
  geom_line(aes(time, value, color = var, linetype = var)) +
  scale_y_continuous(breaks = seq(0, 1, 0.2)) +
  ggsci::scale_color_d3(name = "Stage", labels = Stages) + 
  scale_linetype_discrete(name = "Stage", labels = Stages) +
  ggtitle("Figure 9: A difference model solved with <span style=\"font-family:'Lucida Console', monospace\">method = \"iteration\"</span>",
          subtitle = "Teasel stage distribution") +
  theme(plot.title = ggtext::element_markdown(size = 12),
        panel.grid = element_blank())

10 Plotting deSolve Objects

10.1. Plotting Multiple Scenarios

As an example of plotting multiple outputs on one plot we implement the simple combustion model, which can be found on http://www.scholarpedia.org/article/Stiff_systems: \[ y' = y^2 \cdot (1 - y) \]

The model is run with 4 different values of the initial conditions: \(y\) = 0.01, 0.02, 0.03, 0.04 and written to deSolve objects out1, out2, out3, out4:

combustion <- function (t, y, parms) list(y^2 * (1-y) )

yini <- 0.01
times <- 0 : 200

out1 <- ode(times = times, y = yini, parms = 0, func = combustion)
out2 <- ode(times = times, y = yini*2, parms = 0, func = combustion)
out3 <- ode(times = times, y = yini*3, parms = 0, func = combustion)
out4 <- ode(times = times, y = yini*4, parms = 0, func = combustion)

Plot with ggplot2:

map(list(out1, out2, out3, out4), data.frame) %>% 
  dplyr::bind_rows(.id = "out_") %>% 
  ggplot() +
  geom_line(aes(time, X1, color = out_, linetype = out_)) +
  labs(x = "time", y = "", color = "yini*i", linetype = "yini*i",
       title = "Figure 10: Plotting 4 outputs in one figure") +
  scale_y_continuous(breaks = seq(0, 1, 0.2)) +
  ggsci::scale_color_d3() +
  theme(legend.position = c(1, 0),
        legend.justification = c(1, 0),
        legend.background = element_blank(),
        panel.grid = element_blank())

10.4. Plotting multi-dimensional output

1-D model output

We demonstrate this with a Lotka-Volterra model, implemented in 1-D. The model describes a predator and its prey diffusing on a flat surface and in concentric circles. This is a 1-D model, solved in the cylindrical coordinate system. Note that it is simpler to implement this model in R package ReacTran.

We start by defining the derivative function:

lvmod <- function(time, state, parms, N, rr, ri, dr, dri) {
  with (as.list(parms), {
    PREY <- state[1:N]
    PRED <- state[(N+1):(2*N)]
    
    ## Fluxes due to diffusion
    ## at internal and external boundaries: zero gradient
    FluxPrey <- -Da * diff(c(PREY[1], PREY, PREY[N]))/dri
    FluxPred <- -Da * diff(c(PRED[1], PRED, PRED[N]))/dri
    
    ## Biology: Lotka-Volterra model
    Ingestion    <- rIng * PREY * PRED
    GrowthPrey   <- rGrow * PREY * (1-PREY/cap)
    MortPredator <- rMort * PRED
    
    ## Rate of change = Flux gradient + Biology
    dPREY <- -diff(ri * FluxPrey)/rr/dr + GrowthPrey - Ingestion
    dPRED <- -diff(ri * FluxPred)/rr/dr + Ingestion * assEff - MortPredator
    
    return (list(c(dPREY, dPRED)))
  })
}

Then we define the parameters, which we put in a list:

R <- 20                        # total radius of surface, m
N <- 100                       # 100 concentric circles
dr <- R/N                      # thickness of each layer
r <- seq(dr/2,by = dr,len = N) # distance of center to mid-layer
ri <- seq(0,by = dr,len = N+1) # distance to layer interface
dri <- dr                      # dispersion distances
parms <- c(Da = 0.05,          # m2/d, dispersion coefficient
           rIng = 0.2,         # /day, rate of ingestion
           rGrow = 1.0,        # /day, growth rate of prey
           rMort = 0.2 ,       # /day, mortality rate of pred
           assEff = 0.5,       # -, assimilation efficiency
           cap = 10)           # density, carrying capacity

After defining initial conditions, the model is solved with routine ode.1D:

state <- rep(0, 2 * N)
state[1] <- state[N + 1] <- 10
times <- seq(0, 200, by = 1) # output wanted at these time intervals
print(system.time(
  out <- ode.1D(y = state, times = times, func = lvmod, parms = parms,
                nspec = 2, names = c("PREY", "PRED"),
                N = N, rr = r, ri = ri, dr = dr, dri = dri)
))

   user  system elapsed 
  0.333   0.008   0.340 

The summary method provides summaries for both 1-dimensional state variables:

summary(out)

Extract the data for PREY and PRED:

prey <- subset(out, select = "PREY") %>% 
  data.frame() %>% 
  cbind(time = 0:200) %>% 
  dplyr::select(time, everything())

pred <- subset(out, select = "PRED") %>% 
  data.frame() %>% 
  cbind(time = 0:200) %>% 
  dplyr::select(time, everything())

Plot in 2-d using ggplot2:

# Create a data frame of vars (X1, X2, etc.) and r values for joining to the output:
var_r_prey <- data.frame(
  var = paste0("X", 1:length(r)),
  r
)

var_r_pred <- data.frame(
  var = paste0("X", 101:(101+length(r)-1)),
  r
)

# Pivot the output to longer format and join with the above data frame to add the Distance values:
prey_long <- data.frame(prey) %>% 
  pivot_longer(cols = starts_with("X"), names_to = "var") %>% 
  left_join(var_r_prey, by = "var")

pred_long <- data.frame(pred) %>% 
  pivot_longer(cols = starts_with("X"), names_to = "var") %>% 
  left_join(var_r_pred, by = "var")

# Plot:
P14a <- ggplot(prey_long) +
  geom_raster(aes(time, r, fill = value),
             interpolate = TRUE) +
  geom_contour(aes(time, r, z = value), bins = 12,
               color = "white", size = 0.1, alpha = 0.5) +
  scale_fill_distiller(palette = "Spectral",
                       breaks = seq(0, 10, 1)) +
  labs(x = "time", y = "", title = "PREY", fill = NULL) +
  guides(fill = guide_colourbar(barheight = unit(2.7, "in"),
                                ticks.colour = "black",
                                frame.colour = "black")) +
  coord_cartesian(expand = F) +
  theme(plot.title = element_text(hjust = 0.5))

P14b <- ggplot(pred_long) +
  geom_raster(aes(time, r, fill = value),
             interpolate = TRUE) +
  geom_contour(aes(time, r, z = value), bins = 12,
               color = "white", size = 0.1, alpha = 0.5) +
  scale_fill_distiller(palette = "Spectral",
                       breaks = seq(0, 10, 1)) +
  labs(x = "time", y = "", title = "PRED", fill = NULL) +
  guides(fill = guide_colourbar(barheight = unit(2.7, "in"),
                                ticks.colour = "black",
                                frame.colour = "black")) +
  coord_cartesian(expand = F) +
  theme(plot.title = element_text(hjust = 0.5))

P14a + P14b + plot_annotation(title = "Figure 14: 2-d plots of PREY and PRED") + plot_layout(guides = "collect")

Plot PREY in 3-d using plotly:

# Get the output into a list suitable for plotting a surface plt in plotly:
prey_3d <-
  list(
    x = times,
    y = r,
    z = t(subset(out, select = "PREY"))) # need to transpose this to match the x as times and y as r

# ColorBrewer Spectral palette:
my_palette <- rev(colorRampPalette(RColorBrewer::brewer.pal(9, "Spectral"))(40))

# Plot:
plot_ly(x = prey_3d$x, y = prey_3d$y, z = prey_3d$z, colors = my_palette) %>%
  add_surface() %>%
  layout(
    scene = list(
      camera=list(
        eye = list(x=0, y=-2, z=1)
        )
    )
  )

Plot PRED in 3-d using plotly:

pred_3d <-
  list(
    x = times,
    y = r,
    z = t(subset(out, select = "PRED"))) # need to transpose this to match the x as times and y as r

# Plot:
plot_ly(x = pred_3d$x, y = pred_3d$y, z = pred_3d$z, colors = my_palette) %>%
  add_surface() %>%
  layout(
    scene = list(
      camera=list(
        eye = list(x=0, y=-2, z=1)
        )
    )
  )

2-D model output

Consider the very simple 2-D model (100*100), containing just 1st order consumption, at a rate r_x2y2, where r_x2y2 depends on the position along the grid. First the derivative function is defined:

Simple2D <- function(t, Y, par) {
  y <- matrix(nrow = nx, ncol = ny, data = Y) # vector to 2-D matrix
  dY <- - r_x2y2 * y # consumption
  
  return(list(dY))
}

Then the grid is created, and the consumption rate made a function of grid position (outer):

dy <- dx <- 1 # grid size
nx <- ny <- 100
x <- seq (dx/2, by = dx, len = nx)
y <- seq (dy/2, by = dy, len = ny)

# in each grid cell: consumption depending on position
r_x2y2 <- outer(x, y, FUN=function(x,y) ((x-50)^2 + (y-50)^2)*1e-4)

After defining the initial values, the model is solved using solver ode.2D. We use Runge-Kutta method ode45:

C <- matrix(nrow = nx, ncol = ny, 1)
ODE3 <- ode.2D(y = C, times = 1:100, func = Simple2D, parms = NULL, dimens = c(nx, ny), names = "C", method = "ode45")

Print a summary::

summary(ODE3)

Plot using ggplot2:

expand_grid(x, y) %>% 
  mutate(z = as.numeric(r_x2y2)) %>% 
  ggplot(aes(x, y, z = z)) +
  geom_text_contour() + # requires the metR package
  geom_contour() +
  scale_x_continuous(breaks = seq(0, 100, 20)) +
  scale_y_continuous(breaks = seq(0, 100, 20)) +
  ggtitle("Figure 15: Contour plot of 2-D variables") +
  coord_fixed()

---

1. Soetaert K, Herman PMJ (2009). A Practical Guide to Ecological Modelling. Using R as a Simulation Platform. Springer. ISBN 978-1-4020-8623-6.↩︎

2. Shampine L, Thompson S (2000). Solving Delay Differential Equations with dde23. URL http://www.runet.edu/~thompson/webddes/tutorial.pdf.↩︎